assign non-target gopher-black-and-yellow rockfish with rubias

To start off with, let’s load data and libs:

Some initial filters

Take highest read-depth call for multiply-genotyped DNA_IDs

I’m not sure if there are any of these, but best to leave it in here… if an individual is multiply-genotyped, take the genotype with the highest total read depth.

Remove the 6 loci which Hayley has been removing

# read in a list of the 6 loci
to_remove <- read_csv("../data/loci_to_remove.csv")

Parsed with column specification:
cols(
  locus = col_character()
)

# only keep the loci that are not those 6
keepers <- geno_one_each %>%
  anti_join(., to_remove, by = "locus")

that should leave 90 loci

Toss out indivs with data at fewer than 25 loci

Now, toss out any individual with fewer than 25 non-missing loci

no_hi_missers <- keepers %>% 
  group_by(NMFS_DNA_ID) %>%
  filter(sum(!is.na(allele)) >= (25*2))

So, we started with 12852 and after filtering out indivs with fewer than 75 genotyped loci, we were left with 12741 individuals. Those are the ones that we will run through rubias to identify to species.

Read in baseline genotypes and remove loci and individuals with too much missing data

From the .csv output from microhaplot, we need to reformat the data for rubias.

# add reference column to prepare data for rubias
spp.id_loc1 <- spp.id_loc %>%
  mutate(sample_type = "reference")
# add the repunit column that is identical to the group in this case
x <- spp.id_loc1 %>%
  mutate(repunit = group)
# reorder the columns and get the data frame in the right format
spp.id1 <- x[,c(6,7,1,3,2,4:5)]
spp.id2 <- spp.id1 %>%
  rename(collection = group) %>%
  rename(indiv = indiv.ID)

To ensure that haplotypes are coded identically across both the reference and mixture genotypes, we need to combine the two data sets, and then change haplotypes into numeric alleles.

str(alleles)

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   2548448 obs. of  7 variables:
 $ sample_type: chr  "reference" "reference" "reference" "reference" ...
 $ repunit    : chr  "Sauriculatus" "Sauriculatus" "Sauriculatus" "Sauriculatus" ...
 $ collection : chr  "Sauriculatus" "Sauriculatus" "Sauriculatus" "Sauriculatus" ...
 $ indiv      : chr  "s1" "s1" "s1" "s1" ...
 $ locus      : chr  "Plate_1_A01_Sat_GW603857_consensus" "Plate_1_A11_Sat_GE820299_consensus" "Plate_2_A09_Sat_EW986980_consensus" "Plate_2_C08_Sat_EW987116_consensus" ...
 $ gene_copy  : num  1 1 1 1 1 1 1 1 1 1 ...
 $ allele     : chr  "ACTAGCGTT" "CGGCTCCGGGTATCT" "CGAGATAGXGTCTGTT" "AGACTTCGG" ...

Here we turn alleles into integers, spread the data frame, and then get it into the right format to run rubias:

We will use infer_mixture in rubias, which requires two separate data frames, one with the reference genotypes and the other with the mixture.

I’ll split the data frame that I created (but it needed to be bunged together for the conversion of alleles to integers).

Now, with the two data frames separated, run the infer mixture analysis:

The posterior means of group membership in each collection is in the PofZ column.

Make a quick little plot to display the diversity of species in our NSF juvenile samples.

# ordered by number of samples
ggplot(grp_tally2, (aes(x = reorder(repunit, -n), y = n))) + 
  geom_bar(stat = "identity") +
  theme_bw() +
  ylab("Samples") +
  xlab("Species") +
  theme(
    axis.text.x = element_text(angle = 45, size = 16),
    axis.text.y = element_text(size = 16),
    axis.title = element_text(size = 18)) +
  scale_color_brewer(palette = "Set1")
  
ggsave("pdf_outputs/nsf_juv_id.pdf", height = 7, width = 12)

# There are individual assignments for each sample to each reference.
# If I want to keep the top assignment for each sample:
kept_assignments <- mix_assign$indiv_posteriors %>%
  group_by(indiv) %>%
  filter(PofZ > 0.5) %>%
  arrange(desc(PofZ))
kept_assignments

Picking out the gopher, black-and-yellow, and copper rockfishes

Interesting to see all the other species in there - mostly china and quillback. (And in reality, these non-target species only make up 1-2% of samples.)

But what we actually want is a list of NMFS IDs for the fish assigned to gopher, black-and-yellow, and copper, which we can then use to select those fish from our NSF dataset, where we presumably have many more alleles for those species because the VCF used to call genotypes included hundreds of samples for each of those species.

This dataset includes 12,741 samples.

That’s 7,519 samples.